Dynamical Instabilities of Quasi-static Crack Propagation Under Thermal Stress 
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We address the theory of quasi-static crack propagation in a strip of glass that is pulled from a hot 
oven towards a cold bath. This problem had been carefully studied in a number of experiments that 
offer a wealth of data to challenge the theory. We improve upon previous theoretical treatments in a 
number of ways. First, we offer a technical improvement of the discussion of the instability towards 
the creation of a straight crack. This improvement consists of employing Pade approximants to 
solve the relevant Weiner-Hopf factorization problem that is associated with this transition. Next 
we improve the discussion of the onset of oscillatory instability towards an undulating crack. We offer 
a novel way of considering the problem as a sum of solutions of a finite strip without a crack and an 
infinite medium with a crack. This allows us to present a closed form solution of the stress intensity 
factors in the vicinity of the oscillatory instability. Most importantly we develop a dynamical 
description of the actual trajectory in the regime of oscillatory crack. This theory is based on the 
dynamical law for crack propagation proposed by Hodgdon and Sethna. We show that this dynamical 
law results in a solution of the actual track trajectory in post critical conditions; we can compute 
from first principles the critical value of the control parameters, and the characteristics of the 
solution like the wavelength of the oscillations. We present detailed comparison with experimental 
measurements without any free parameter. The comparison appears quite excellent. Lastly we show 
that the dynamical law can be translated to an equation for the amplitude of the oscillatory crack; 
this equation predicts correctly the scaling exponents observed in experiments. 
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I. INTRODUCTION 

In 1993 Yuse and Sano reported a simple experiment 
on fracture in glass that nevertheless has attracted 
great attention from the fracture community. The experi- 
ment examined a strip of glass pulled at constant velocity 
v from an oven into water, cf. Fig. ^ At small enough 
velocity nothing happens. A first critical velocity heralds 
the onset of a straight crack, whereas exceeding a second 
critical velocity results in an oscillatory crack. Finally, 
at sufficiently high velocities the crack pattern exhibits 
multiple fractures and disorder. The reason for the high 
interest in this relatively simple experiment is of course 
that it offers a challenge for the theoretical description 
of fracture processes. Being essentially a "quasi-static" 
process, as the velocity v is very much smaller than the 
Rayleigh speed, the fracture process here is free of many 
of the complications arising in truly dynamic fracture |2| • 
Nevertheless, in the absence of a microscopic theory of 
the "process zone" (how materials actually break) even 
the dynamics of quasi-static crack propagation in brittle 
materials remains a debatable issue. 

The lack of dynamical theory for the cracking pro- 
cess does not hinder the understanding of the onset of 
straight cracks in the above experiment. Indeed, already 
in the year following the original experimental observa- 
tions, Marder has set up the equations describing the ef- 
fect of the temperature field on the elastic theory of the 
material, and presented a qualitative description of the 
onset of straight cracks Q- From the quantitative point 
of view this treatment was lacking, in particular the frac- 
ture energy turned out to be strongly velocity dependent 
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FIG. 1: Schematic representation of the experiment: a thin 
glass plate is pulled at a velocity v away from a heater into 
a cold bath. The control parameters are the temperature 
difference between the oven and the water AT, the pulling 
velocity v, the spatial separation between the thermal baths 
h and the width of the plate 2b. 



against physical intuition. The tools employed by Marder 
did not allow however a correct prediction of the charac- 
teristics of the oscillatory crack propagation. The next 
decisive theoretical step was taken by Adda-Bedia and 
Pomeau Q. These authors reproduced Marder's results 
for straight cracks, but also developed a successful cri- 
terion for the secondary instability to oscillatory cracks. 
They employed the universal form of the near-tip stress 
tensor field, i.e. 



<7y(r,l 



(1) 



Here Ki and Ku are the "stress intensity factors" with 
respect to the opening and shear modes, whereas 
and Yifj(Q) are universal angular functions common to 
all configurations and loading conditions. Adda-Bedia 
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and Pomeau invoked the well-known and extensively used 
"principle of local symmetry" , which states that the path 
taken by a crack in brittle homogeneous isotropic mate- 
rial is such that the local stress field at the tip of the crack 
is of mode I type, annulling Kjj. The considerations in 
led to the conclusion that the appearance of a negative 
Kji for positive deviations from straight trajectory (or 
positive Kn for negative deviations) was tantamount to 
the onset of the oscillatory instability. Nevertheless these 
authors did not offer a careful quantitative comparison 
against the experiments knwon at the time. Their predic- 
tion of the fracture energy and the wavelength of oscilla- 
tions differed significantly from the experimental values. 

In light of these results, we should explain at this point 
why do we feel that further theory is called for. First, we 
point out that "principle of local symmetry" is hardly a 
dynamical theory. It can predict an instability, but taken 
literally would only agree with a fracture path that has 
sharp kinks. It cannot be employed to predict the actual 
trajectory of a slowly moving crack when the latter is not 
straight. Second, since the theoretical works cited above 
there have been additional experimental studies of this 
very same problem [a, [J la] , offering a wealth of data to 
challenge the theory, a challenge that had not been picked 
up by the theorists. Last, but not least, we feel that we 
can improve on a number of technical issues tackled by 
previous authors; these will be spelled out in the sequel, 
hopefully gratifying the diligent reader as we go along. 

From the conceptual point of view we offer a point of 
departure from previous treatments by adopting a dy- 
namical description of the crack development. In this we 
follow Hodgdon and Sethna who have built upon the 
principle of local symmetry, using standard theoretical 
methods, to reach a dynamical law for crack propagation 
which is given by 



^ = -/ K U n (2) 



recent mathematical advances |10| we use Pade approxi- 
mants to significantly improve the treatment. Since the 
results of this improved treatment are relied upon in our 
solution of the oscillatory instability, we present the the- 
ory for the straight crack in some detail. Sect. ^Q] intro- 
duces the main novel results of our study in the context of 
the secondary instability to oscillatory cracks. Using the 
dynamical law Eq. Q we show that as one crosses the 
second critical value of the parameters the solution of the 
equations changes its nature. We solve the equation near 
the onset of the oscillatory instability and calculate the 
critical values of the control parameters, the wavelength 
of oscillations and the material function /. Although 
our criterion for the oscillatory instability is in agree- 
ment with 

HE1 (which was based on the "principle of 
local symmetry"), we can go considerably further in de- 
scribing the actual dynamics in the oscillatory regime. In 
particular we present a quantitative comparison with the 
experiments. Our handling of the oscillatory instability 
includes also a technical improvement on the analysis of 
; the latter needed a separate Weiner-Hopf problem for 
every order in the amplitude of the perturbation. In our 
calculation we derive a new expression for Kn to leading 
order in the amplitude of the oscillations, an expression 
that requires a solution of only one Wiener-Hopf problem. 
This simplification is achieved by presenting a new way 
to decompose the straight crack problem into a singular 
and non-singular parts and then using a classical result 
of Cotterell and Rice . A crucial step in the calcula- 
tion is the Wiener-Hopf factorization, for which we apply 
the new method of solution based on Pade approximants 
|ldj | . Employing the dynamical law of crack- tip propaga- 
tion we calculate the critical exponents for the transition 
and compare them with the experimental data. Sect. IIVI 
offers a summary and concluding remarks. 



II. THE STRAIGHT CRACK 
A. Preliminaries 



where i and n are the tangent and the normal to the 
crack tip respectively, and / > is a material parameter 
that we assume to be nearly independent of t and n in 
the quasi-static limit. This law predicts a differ entiable 
crack path such that Kjj is reduced. We will demon- 
strate that this law of motion provides us with predictions 
which are in excellent agreement with the characteristics 
of the crack trajectory in the oscillatory regime. We be- 
lieve that this is the first context in which Eqs. (J2J are 
compared against a challenging set of experimental data; 
The comparison appears quite favorable. 

In Sect. [H]we discuss, facing the danger of being super- 
fluous, the problem of the primary instability leading to 
a starting crack propagation once more. This instability 
had been correctly treated in 0, 3 , but we offer a techni- 
cal improvement in the handling of the Weiner-Hopf fac- 
torization problem that is the basis of the solution. Using 



By varying the experimental control parameters one 
varies the amount of elastic energy stored in the glass 
plate. One can choose various paths in parameter space; 
in this work we adopt the scheme of Q , fixing the values 
of AT, h and v. The growth state depends then on the 
plate's width 2b: for small enough values of the width a 
seeded crack does not grow; For a width greater than a 
critical value L c , a crack, whose tip penetrates a length 
£ away from the cold bath, moves at a velocity —v. This 
crack is stationary in the laboratory frame of reference 
and is stable as long as the width is smaller than another 
critical value L osc . Above this value the crack becomes 
unstable and exhibits an oscillatory lateral motion with 
a well-defined amplitude and wavelength, still traveling 
at a velocity —v. As the width is further increased the 
propagation becomes less and less regular. 

The no crack - straight crack propagation transition is 
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well-understood and the agreement with the experimen- 
tal data is favorable 6]. In this case the propagation is 
pure mode I and the transition is governed by the follow- 
ing Irwin's relation 



Kl 
E 



(3) 



where V is the fracture energy which is a material prop- 
erty and E is Young modulus. In this section we address 
again this problem and cite the formal solution. The 
merit of our treatment will be in providing a detailed 
scheme for performing the Wiener-Hopf factorization in 
a novel way. 



B. The formulation of the problem 

Imagine the glass plate as in Fig. ^ with a straight 
crack penetrating into the glass from the water side. We 
choose a coordinate system such that x — is at the 
crack tip, (marking the water level at x = —£ where 
£ is the penetration depth of the straight crack). The 
y coordinate spans the interval [—6,6]. The condition 
for mechanical equilibrium under plane strain conditions 
caused by a nonuniform temperature field reads 



V 2 V 2 X (x,y) = -Ea T V 2 T(x) , 



(4) 



where ax is the thermal expansion coefficient and x IS 
the Airy potential which is related to the stress tensor 
by 



dy 2 ' 



dx 2 ' 



d 2 x 
dxdy 



Using the symmetry of the problem we state the bound- 
ary conditions as follows: 

<7 yy (x, ±6) = a xy ( x , ±6) = o X y (x , 0) = (6) 



Uyy (x, 0) = for x < 0, u v (x, 0) = for x > . (7) 

Fourier transforming Eq. 10} in the x direction and focus- 
ing on the upper half plate one obtains 3] the following 
Wiener-Hopf equation |12| 



with 



. (k, 0) = -F (k) u y (k, 0) + D e (k) 



F{k) = Ek s ™ h2 ( kb )~ kb2 



(8) 



sinh(2kb) + 2kb 



(9) 



D e (k) = 2Ea T fi (fc) 



[1 - cosh(kb)] [sinh(kb) - kb] 



sinh(2kb) + 2kb 



where one still has to obey the boundary conditions of Eq. 
. Note that the subscript I denotes the transformation 
x — > x + £ in the temperature field such that the origin 
of the coordinates system is at the tip of the crack. For 
convenience, from now on, we rescale all lengths in the 
problem by the half width b. 

Writing F (k) = , where F~ (k) has neither zeros 

nor singularities for Im (fc) < 0, and F + (k) has neither 
zeros nor singularities for Im (k) > 0, the Wiener-Hopf 
method [Ta| results in 



Uy (fc,0): 



1 



F- (k). 





dx 

-oo 



dk 

2^ 



D e (k)F+(k)e 
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(11) 



From this solution one can extract [4| the mode I stress 
intensity factor introduced in Eq. ijljl. 

f°° dk 

K I (AT,v,b 1 h)= De (k)F+(k) (12) 

J-oo 27T 

Note that if the fracture energy T is known then the mode 
I stress intensity factor characterizes completely the no 
crack - straight crack transition. Kj is a positive quantity 
and is different from zero, as a function of £, only on a 
scale of the order of b. To calculate this quantity we 
need first to solve for the temperature field and second 
to provide a method for accomplishing the Wiener-Hopf 
factorization. 



C. The temperature field 

The nonuniform temperature field induces the stress 
field in the elastic plate. In this subsection we solve for 
the temperature and study its properties. For simplicity, 
we set the zero of the coordinates system at the cooling 
front level to avoid £ dependence which is unnecessary 
in the present context. For later calculations we will use 
the aforementioned transformation to put back the £ de- 
pendence. 

In the frame of reference of the plate the temperature 
field obeys the heat equation 



^ = DV 2 T 
dt 



with the boundary conditions 

VT-n = 
T(x = 0) = T 
T(x = h) = T + AT. 



(13) 



(14) 



(10) 



Here D ~ 0.47mm 2 / sec is the diffusion coefficient of the 
glass, h is the distance between the cold bath and the 
heater and n is the unit vector normal to the boundary 
of the domain. 
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between the thermal baths, but assumed perfect thermal 
baths, an assumption that leads to a discontinuity of the 
gradient near the baths. This discontinuity results in 
incorrect estimates of the stress field; since at low veloci- 
ties (see Fig. 2) the only gradient variations are in these 
regions, we cannot expect a very good quantitative agree- 
ment with the experiment for low velocities. At higher 
velocities, there is a significant curvature inside the glass, 
so we expect a better quantitative agreement with the ex- 
periment. In ||| the temperature field was measured and 
found to vary smoothly near the baths due to the finite 
impedances of the baths. In the absence of the exper- 
imental data of the measured temperature field we will 
use the ideal baths approximation. 



position (mm) 



FIG. 2: The calculated temperature distribution inside the 
glass plate. We can identify two distinct regimes, the diffusive 
regime at low velocities (v = 0.05mm/sec in the lower curve), 
and the advective regime at higher velocities (v = 0.3mm/sec, 
and v — 0.5mm/ sec in the upper two curves). These curves 
should be compared with the measured temperature field in 



This equation can be simplified for the straight crack 
configuration, for which there is no y dependence, by 
looking for a stationary solution in the laboratory frame 
of reference, of the form T(x — vt). This solution obeys 
the stationary diffusion equation 



D. The Wiener- Hopf factorization 

The crucial element in the solution is the Wiener-Hopf 
factorization of the known kernel F (k). Generally 
it is not possible to find an exact factorization, so 
one tries an approximant F(k) ~ F(k) that can be 
exactly factorized. This can be made rigorous following 
a recently proven theorem jlfj which establishes the 
closeness of the product factors of F (k) to those of 
F (k) in their region of regularity if F (k) ~ F (k) for 
all k G D, where D is the strip of analyticity of F(k). 
A commonly followed first step in finding F(k) is to 
examine the behavior of F(k) near zero and ±oo 



V 2 T 



1 dT 

d t h dx 



= 



(15) 



where dth = % is the thermal diffusion length. The 
exact solution of this equation is 



k 

F(k) — > ± — as k — > ±oo 



F(k) 



k± 

12 



as k — * ±0 



(17) 
(18) 



T (x) = AT 



-x/dt 



1 



-h/dt. 



'(x)9(h-x)+0(x- h) 



(16) 

We can identify two distinct regimes, as shown in Fig. 
2. The first one, at low velocities, is the diffusive regime, 
in which the temperature field is controlled by the spa- 
tial separation of the thermal baths h. The second one, 
at higher velocities, is the advective regime, in which 
the temperature field is controlled by the thermal dif- 
fusion length dth- Actually, there is a third regime, at 
still higher velocities, in which the temperature field is 
controlled by the thickness of the plate. Note that we as- 
sumed that the temperature is uniform along this dimen- 
sion and therefore the breakdown of this assumption in 
this regime leads to a three-dimensional problem, which 
is outside the scope of our two-dimensional model. 

The temperature field enters the problem through the 
term —Eolt^ 2 T (x) of the inhomogeneous Bi-Laplace 
equation. This term is sensitive only to variations of the 
temperature gradient, i.e. to the curvature of the thermal 
field. In Eq. i|ltj|) we considered a finite spatial separation 



A standard approach to finding a good factorization is 
then to seek a function <fi(k) that reproduces the asymp- 
totic behavior of F(k) and to correct it by a ratio of two 
polynomials 



F(k) = 
where for example 

<f>(k) 



fc 4 + ak 2 + p 
k 4 + -fk 2 +p 



\/4F + 144 



(19) 



(20) 



and a, j3, 7 are free parameters that should be chosen 
as to best fit F (k) to F(k). In principle, one can use 
higher order polynomial ratio to achieve greater accuracy. 
The disadvantage of this approach is that the positions 
of the poles and zeros are not well-controlled and that 
the convergence behavior of the process is not clear. 

In our work we follow a new method developed in 
[lOj. In the heart of this approach lies the use of Pade 
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approximants. A approximant of F (k) is written 

as 



F(k) 



P N {k) 
Qm (k) ' 



(21) 



where 



P N (k) = a + aifc + a 2 k 2 H ^ a N k N (22) 

(fe) = l + 6 1 fc + 6 2 fc 2 + --- + 6 M fc M (23) 



The coefficients (2yj . ^yj CITS determined from the Taylor- 
series expansion of F (k) at any regular point. Let us take 
the expansion point to be k = 0, so 



F(k) = ^c„fc" 



(24) 



71=0 



where c„ are known. In order to solve for the unknown 
coefficients one should set 



P N (k) 
(k) 



+ 0{k N+M+1 )=J2cnk r > 



(25) 



n=0 



to obtain a set of linear equations 01 ■ I n this method 
one approximates directly the factorization and the 
process is completely algorithmic. Note that since 
in practice one uses the truncated series of F (k) it 
is not possible to approximate directly the product 
factors for arbitrary large k's. In order to overcome 
this difficulty we should find the asymptotic form of 
the factorization and use it as a boundary condition for 
the Pade approximants. The asymptotic factorization 
is found by noticing that the zeros and poles of F(k), 
which are respectively the solutions of the equations 
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FIG. 3: The relative error. 
k. 



F~(fc)/F+(fc)-F(fc) 
F(fc) 



as a function of 



Using this relation and the asymptotic relation of Eq. 
(|17|l we conclude that the asymptotic factorization is 



F + (k) 




(30) 



We should choose the Pade approximants to match 
these asymptotic forms. This is achieved by squaring 
the original kernel F(k), to obtain an even function of k 
that behaves els k cis |fc| — > oo. Hence, we can derive an 
[nv^] P a de approximation 



F 2 (k) 



P 



N+2 



(k) 



>N 



(AO 



(31) 



with N even. This approximation contains N + 2 zeros 
and N poles which become, after taking the square-root, 
N — 1 branch points in the upper half-plane and N — 1 
branch points in the lower half-plane. The relative error, 
r(fc)/ ff" F(t,) , for N = 28 is shown in Fig. 3. The 

asymptotic matching of F + (k) to J is shown in Fig. 



sinh 2 (w„) — w 2 n — 
sinh(2z„) — 2z n = 



(26) 



(27) 



have the property that if, for example, w n is a zero of 
F(k) then w n , —w n and —w n are also zeros. The same 
holds for the poles. Therefore, considering the solutions 
of Eas. (|26l27ll only in the first quadrant, we obtain 



F(k) 



fc 4 aT=i(l-^)(l + ^r)(l 



-)(! + £) 



12 ir=i(i 



From here it follows that 

k 2 n£=i(i 



(28) 



F-{k) 



V12 n~i(l -£)(! + £) F+(-k) 



(29) 



E. The determination of the fracture energy V 



In Sect. IIIII we develop a theory for the oscillatory in- 
stability and compare it with the experimental results. A 
crucial parameter in that theory, as here, is the fracture 
energy T. We can extract the fracture energy from the 
experimental threshold for propagation, cf Eq. © • The 
relevant measurement was reported in |6j in which a lin- 
ear elastic model identical to ours was used to extract the 
fracture energy. T was chosen to best fit the experimen- 
tal data of the onset of the straight crack propagation. It 
was found that T depends weakly on the velocity when 
the idealized thermal profile was employed; T turned out 
velocity independent for the actual thermal profile mea- 
sured in the experiment. Since we do not have the experi- 
mental data for the thermal profile we will use, for consis- 
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FIG. 4: The real and imaginary parts of F + (k) compared to 
the real and imaginary parts of yj -^jj:, as a function of k. 

tency, a typical value of the former, i.e. T ~ 3.75 J/m 2 . 
This value should be compared to Fig. 5 in Ref. J6j. 

III. THE STRAIGHT TO OSCILLATORY 
CRACK TRANSITION 

Solving a dynamic fracture problem in the quasi-static 
limit consists of (i) solving the equilibrium equations for 
the stress field together with a given set of boundary 
conditions at the sample boundaries and on the (a-priori 
unknown and evolving) crack boundary and (ii) employ- 
ing a dynamical principle to evolve the crack. A proper 
solution determines the correct shape of the crack as a 
function of time. Clearly, the predictions of the employed 
crack growth law should be consistent with the experi- 
mental observations. In this section we use the dynamical 
law in Eq. to study the straight to oscillatory crack 
transition. 

Rewriting the tangential and normal unit vectors at 
the tip of the crack in terms of the angle 8 that the tan- 
gential unit vector makes with the a;-axis we obtain 

t = cos 8 x + sin 8 y 
h = — sin 8 x + cos 8 y (32) 

which, upon substitution into Eq. J5J), leads to 

This equation predicts that as long as Kn — 0, the 
crack will propagate in a straight line. Nevertheless, in 
any real material there exist intrinsic instabilities, due to 
imperfections of the material and the loading conditions, 
which produce a small random Kn ^ 0. We are now 



facing two distinct questions: Under what conditions, in 
terms of the control parameters AT, h, v and 2b, the 
straight crack propagation becomes unstable? Once the 
straight crack propagation becomes unstable, what is the 
stable stationary path that it follows? 

The criterion of stability arises naturally from the dy- 
namical equation. If 8 and Kxi have the same sign, with 
/ > 0, then || has the opposite sign and \8\ decreases, 
which means that a small perturbation decays. By the 
same argument, for 8 and Kjj having the opposite sign 
a small perturbation grows. This criterion is identical to 
the one suggested in Refs. 0,0- 

The question of the future evolution of the crack, once 
the instability threshold was reached, should be answered 
by solving the dynamical equation just above the onset of 
the instability. Guided by the experimental observation 
that the shape of the crack just above the onset of the 
oscillatory instability is a pure sine function, we introduce 
a smooth deviation from straight crack path 

y(x,t) ~ Asin[w(x + vt)} + 0(A 3 ) for x < (34) 

This assumption serves two roles: first, it represents a 
single mode component, corresponding to a wavenumber 
W, in the linear decomposition of a small random pertur- 
bation on top of the straight crack and will enable us, for 
t = 0, to analyze its stability; Second, it is an ansatz for 
the solution of the dynamical equation just above the on- 
set of instability. Note that y(x,t) satisfies \y(x, t)\ <C 1 
and \y'(x,t)\ < 1. 

In this section we will study the stability of the straight 
crack as well as the time evolution of the crack after the 
onset of instability by analyzing the dynamical model 
of Eq. (J23). The stability is studied by applying the 
stability criterion derived above for which an expression 
of Kji to leading order in the perturbation amplitude 
in needed. We derive this expression by introducing an 
auxiliary problem, the so-called decomposition problem, 
whose effective solution enables one to significantly sim- 
plify the derivation. The critical point is calculated by 
solving the set of equations governing the transition. The 
time evolution of the crack just above the critical point 
is studied by directly solving the equations of motion for 
the tip of the crack. 

A. The decomposition problem 

In order to study the stability of the straight crack to 
small perturbations we arbitrarily choose t = in Eq. 
(|34|l . This choice will be shown later to be legitimate. 
We are interested first in finding an expression for Ku 
to leading order in the amplitude of the perturbation. We 
begin by formulating an auxiliary problem. The presence 
of a crack, which is usually modeled as a mathematical 
branch cut, introduces the famous square-root singularity 
of the stress field near the tip of the crack. We want 
to represent our problem as a sum of two parts. The 
first part contains no singularity, implying that it is crack 
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free, but it includes the geometry of the problem and the 
thermal field. The second part contains the singularity, 
implying that there is a crack with a given load on it, but 
the domain is infinite. The non-singular part is chosen 
such that it reproduces the required boundary conditions 
on the plate's edges and on the crack. 

Once we obtain the load on the semi-infinite straight 
crack in an infinite medium we can apply the classical 
result of Cotterell and Rice for slightly curved cracks [Tl). 
The mathematical formulation of this decomposition 
problem leads to a set of integro-differential equations 
whose complexity may cast doubt on the usefulness of 
the whole procedure In what follows we will show 
how to avoid these mathematical difficulties and effec- 
tively solve the problem. 

To see that the solution is almost at hand, suppose 
for a moment that we succeeded to solve the problem 
in this way for a straight crack and for a given set of 
the control parameters. The load on a straight crack in 
an infinite domain , which is a fictitious quantity, must 
be a pure mode I load by symmetry. We denote it as 
o~ yy (x, y = 0; £), where we marked explicitly the paramet- 
ric dependence on I. The mode I stress intensity factor 
is given by [l3| 



dxo-lJx,y = 0;£) 



(35) 



completely different mathematical procedures, should be 
identical functions of £. The £ dependence of the second 
expression is given by the phase factor e~ M which imme- 
diately implies that & yy (k, y = 0;£) = <J yy (k, y = 0)e~ tke . 
We conclude, by the uniqueness of the Fourier transform, 
that 



/(*,tf = 0) 



D (k)F+(k) 



D (k)F+(k) 



2 fO 

7T J — C 



dxe 



which effectively solves the auxiliary problem. 

Hence, we have shown how one can use the Wiener- 
Hopf solution for a traction-free straight crack in a fi- 
nite configuration in order to find the effective load on 
a straight crack in an infinite configuration via the so- 
lution of the decomposition problem. We reiterate that 
this load is a fictitious tension on a crack in an infinite 
domain corresponding to a traction-free situation in a fi- 
nite domain. The solution of this auxiliary problem will 
enable us later on to use the powerful tool of the complex 
potential method that is most suitable for an infinite do- 
main problems. Other theoretical treatments that were 
unable to solve this problem led to incorrect predictions 

M. 



Introduce now the x- Fourier transform of<7 yy (x, y — 0; I), 
denoted as cr l vv (k, y = 0;£). With this object in mind we 
rewrite Eq. (|35jl as 



Kj(£) 



^ yy (k,y = 0;£)[ 



dxe 



-ikx 



( 36 ) 

On the other hand, we have calculated the same quantity 
using the Wiener-Hopf technique (see Eq. H2fr . 

Kj(£)= / D (k)e- M F+(k) (37) 

These two expressions for Kj(£), though derived through 



B. The critical point 

The calculation now is straightforward. Let us se- 
lect a local coordinate system {r, 9} at every point on 
the crack, with r being the distance from the point 
and 9 the angle, starting with 9 = for the tan- 
gent. In such coordinates the normal opening stress 
T n (x, y(x,t)) = agg(x,y(x,t)) and the tangential shear- 
ing stress Tt(x,y(x,t)) = a r e(x,y(x,t)). Using the load 
(jyy we can find these stress components for any small 
deviation y(x,t) from the straight crack. Applying these 
loads to the classical result of Cotterell and Rice [Tl| . 
with t — 0, we obtain the following expressions for K\ 
and Kjj to leading order in A (see Appendix A) 



Kr 



° dxal y (k,y = 0;£)_ + o(A2) 



K n (A,w,t = 0) 



-x 

2 r° dx [A sin{wx) a[ v {k,y = 0;£)} 

TT 



(_a;)3/2 



A w Ki + 0(A 3 



(39) 



This result shows that only Ku is changed to first order 
in the amplitude of the perturbation. The expression for 
Ku has the general form derived in [{|. It is the sum of 
two competitive terms which the authors of refer to 



as a "physical" shear stress which is a destabilizing term 
(the first term), and as a "geometric" shear stress, which 
is a stabilizing term (the second term) . We expect to find 
a range of the control parameters for which for every w ^ 
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FIG. 5: Knj A versus the dimensionless wavenumber w. For 
fixed AT, v and h, the curves from top down show the in- 
creasing of the stored elastic energy via the increasing of the 
width of the plate. It is clear that there is a critical width for 
which Ku(w) — 0. 



the second term dominates the first, leading to Kjj > 0, 
which implies a stable straight crack propagation. Thus, 
our stability criterion states that the transition between 
straight and oscillatory crack propagation occurs when 
there exist A,w ^ such that Kn(A,w,t = 0) = 0. 
One concludes that the straight crack to oscillatory crack 
transition is governed by the following set of equations 



Kj(b,e) + Kh(b,e,w) 

E 

K n (b,£,w) 
dK n (b,l,w) 
dw 



Kj(b,e) 



E 



= r 



(40) 



where we made explicit the dependence on b, I and w. 

The first equation is the Irwin's relation which ex- 
presses the energy balance between the elastic energy flow 
to the tip of the crack and the fracture energy needed to 
create a new crack surface. This fracture energy is a 
parameter of our model; in previous applications this pa- 
rameter was optimized for agreement with experiments 
0. We cannot afford such luxury since we have deter- 
mined already the parameter in Sec. Ill El Therefore 
in our comparison with experiments the theory is truly 
challenged, and the agreement will be shown to be very 
satisfactory. 

The second and third equations express the stability 
threshold. In order to characterize quantitatively this 
transition we adopt the experimental scheme of Ref. @ , 
in which AT and h are kept fixed and for a given v the 
critical width for the the onset of oscillations, T osc , is 
found. Fig. 5 shows a typical situation in which Kjj/A 
is plotted as a function of w for different values of the 
width 2b at constant AT, v and h. It is shown that 
as one increases the width a solution for the equation 
Kji(A, w,t = 0) = appears. 

We solved the above set of equations graphically by 
the following procedure. We fixed AT = 135° C and 
h = hmm and for each velocity we changed b until we 
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FIG. 6: The critical width for oscillations Lose = 2fo osc ver- 
sus the driving velocity v for AT = 135° C and h = 5mm. 
We didn't calculated these quantities for still higher veloci- 
ties since this regime is controlled by three-dimensional effects 
which are outside the scope of our theory. The theoretical val- 
ues are connected by the dashed line that was added as a guide 
for the eye. The experimental values that are represented by 
the unfilled squares were extracted from Fig. 15 in Ref. |8|. 
The deviation of L osc from the measured data is due to its 
high sensitivity to the fine details of the temperature field, 
which we approximated, in the absence of the measured one, 
using the ideal baths assumption. It is clear, as predicted in 
Sec. Ill CI that the agreement with the experiment is much 
better within the advective regime than within the diffusive 
regime. 



converged to the solution. Note that first one has to 
solve the first equation for I. Fig. 6 shows the critical 
width for oscillations T osc = 2& osc as a function of the 
driving velocity v. The experimental data reported in 
Ref [3 have been added for comparison. 



The deviation of T osc from the measured data is due 
to its high sensitivity to the fine details of the tempera- 
ture field, which we approximated, in the absence of the 
measured one, using the ideal baths assumption. Never- 
theless, as predicted in Sec. Ill (Jl the agreement with the 
experiment is much better within the advective regime, 
in which the temperature field is controlled by dth, than 
within the diffusive regime, in which the temperature 
field is controlled by h. Our solution here yields also the 
wavenumber of the unstable mode. In [{| this wavenum- 
ber was identified with the wavenumber of the actual 
trajectory in the post-critical conditions without further 
discussion. We find this unsatisfactory; whether or not 
this wavenumber will be observed in the actual crack tra- 
jectory also in post-critical conditions depends on the dy- 
namics. If the "fastest growing mode" becomes stabilized 
by the nonlinear terms, then this wavenumber would be 
observed. To assess this issue we must turn next to the 
weakly non-linear theory in the post-critical regime. 
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C. The post-critical time evolution of the crack 

Crack propagation laws used in the literature so far 
were unable to predict analytically non-trivial trajecto- 
ries left behind a crack tip. The set of well-controlled 
experiments described here offers a challenge to any dy- 
namical law, especially near the critical point where the 
crack exhibits a lateral oscillatory motion with a well- 
defined wavelength and amplitude. In this subsection 
we will show that the adopted dynamical law meets that 
challenge. The stationary stable path of the crack just 



above the onset of the oscillatory instability is determined 
by the solution of the dynamical equation near the tran- 
sition. Noting that under our assumptions 9 ~ y'(0,t) 
we obtain 

= -fKu(A, w, t) + 0(A 3 ) (41) 

at 

Deriving the time dependent expression for Kn (see Ap- 
pendix A) and substituting our ansatz 1)34(1 ■ Eq. 141|) 

becomes 



Aw 2 v 



sin(-uwi) = \ — 



dx {2Asui[w(x + vt)}al Jx,y = 0;£)}' 



Asm(wvt)a f (x,y = 0;£) 1 
— h - Aw cos(wvt) Kj 

[42] 



This equation has a trivial solution, i.e. A = 0, which is 
the straight crack. An explicit calculation with A ^ de- 
termines that the RHS of this equation is a pure sine func- 
tion. Thus our ansatz (|34[> can be an actual solution only 
if we can choose the control parameters such as to set the 
phase of the sine function to zero at t — (cf. the LHS). 
We see that this is possible with Ku(A,w,t = 0) = 
which is exactly what was calculated above in the con- 
text of linear stability analysis. Thus if this condition 
can be met, (and if A remains small above the critical 
point, cf. the next subsection) we can indeed identify 
the aforementioned wavenumber as the wavenumber of 
the oscillations in the close vicinity of the critical point. 
We conclude that the equation of motion (|53*|l is consis- 
tent with a pure sinusoidal trajectory, which is an exact 
solution of the post-critical dynamics. This is in good 
agreement with the experimental observations. This re- 
sult also shows that the arbitrary choice t = in the 
linear stability analysis is legitimate since one has to fix 
the phase only at one time point. 

Fig. 7 shows the wavelength of the oscillations A osc = 
2Tr/w osc as a function of the driving velocity v. The criti- 
cal width for oscillations L osc — 2b osc , first shown in Fig. 
6, was added for completeness, while the experimental 
data reported in Q have been added for comparison. It 
is clear that the wavelength of oscillations agrees rather 
well with the experimental data, which confirms the as- 
sertion in Ref. 8] that the oscillation wavelength A osc 
seems far less sensitive to the fine details of the temper- 
ature field than L osc . 

As we have indicated before there are several relevant 
length scales in the problem. Here we have calculated a 
new quantity which has the dimension of length and we 
wanted to explore its dependence on the various length 
scales in the problem. Therefore, we have calculated 
the dimensionless oscillation wavelength X osc /L osc at the 
threshold of instability as a function of the dimension- 
less thermal diffusion length dth/L osc . The results are 
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FIG. 7: The wavelength of the oscillations A osc = 27r/w osc 
versus the driving velocity v for AT = 135° C and h = 5mm. 
The theoretical values are connected by the dashed lines that 
were added as guides for the eye. The experimental points 
that are represented by the unfilled points were extracted from 
Fig. 15 in Ref. The results of Fig. 6 were superimposed 
for completeness. The wavelength of oscillations A osc seems 
far less sensitive to the fine details of the temperature field 
than L osc . 



shown in Fig. 8. We have found that within the advec- 
tive regime, which corresponds to relatively high veloc- 
ities, this function can be well fitted by the linear scal- 
ing law A osc /L osc ~ a + (3d th /L osc with a = 0.12 and 
(3 = 2.1, to be compared with the experimental values 
of a = 0.15 and /3 = 2.5 | and the FEM simulation 
values of a — 0.14 and (3 = 2.1 15]. This result shows 
that the scaling between A osc and L osc is controlled by 
the thermal diffusion length. 

Up to now we have dealt with the matching of the 
phases of both sides of Eq. (|42|l . The matching of the 
amplitudes will enable us to calculate the material func- 
tion /. The RHS of this equation, at the critical point, 
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FIG. 8: Dimensionless oscillation wavelength A osc /L osc at the 
threshold of instability as a function of the dimensionless ther- 
mal diffusion length dth/L osc . It is seen that within the ad- 
vective regime the values can be well fitted by a straight line. 
This result shows that the scaling between A osc and L osc is 
controlled by the thermal diffusion length. 



has the form 



AEctT AT \/ b osc A* sin(w Q 



(43) 



where A* is a dimensionless amplitude to be calculated. 
By equating the amplitudes of both sides of Eq. (|4*2*)l we 
obtain 



/ = 



(44) 



The material function / determines the decay length of 
perturbations with finite Ku back to a pure mode I prop- 
agation in the straight crack regime. A typical length is 
constructed from v osc /g, where g = f EotT&T^Jb osc . Fig. 
9 shows v osc /g (in units of 6 osc ) as a function of v osc 
for velocities in the advective regime where the theory 
agrees with the experiment. It is seen that this length 
(in units of 6 osc ) decreases as the velocity increases. As- 
suming that this behavior is not sensitive to the details 
of the temperature field, it is a challenge to any theory 
that will suggest an independent derivation of /. 



D. The critical exponents of the amplitude of 
oscillations 



In the previous subsections we discussed how the so- 
lution of the equation of motion changes its nature from 
an A = solution to an A solution, but since we 
considered all the relevant quantities to O(A) we could 
not study the properties of the amplitude itself. In or- 
der to extend our analysis we should introduce a time 
dependent amplitude in our ansatz 
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FIG. 9: v osc /(g b osc ) as a function of iw for velocities in the 
advective regime, where g = fEaTAT\/b OEC . 



and write the equation of motion to 0(A 3 ). The orien- 
tation of the tip of the crack is given by 



8{t)=tg- x [y'{x = 0,t)]~y'(x = 0,t) 



y' 3 {x = 0,t) 



which upon substitution of our ansatz becomes 



(46) 



6{t) ~ w A(vt) cos(wvt) + A'(vt) sm(wvt) (47) 
— —[w A(vt) cos(wvt) + A'(vt) sm(wvt)} 3 

o 

where the prime denotes a derivative with respect to the 
argument vt. Because of the symmetry A — > —A the 
next order term in the expansion of Ku in powers of A 
is of 0(A 3 ). Thus, 



K 



ii 



K${A{x + vt) sm[w(x + vt)}} 



-(3) 



(48) 



+ K u '{A(x + vt) sin[w(x + vt)}} + 0{A b ) , 



where -Kr^{-} is the functional that was calculated in the 

frt\ 

Appendix, being of 0(A). K}/{-} is a functional whose 
derivation is straightforward but very lengthy; we do not 
present it explicitly here, but note that it yields a term 
ofO(A 3 ). 

In order to proceed we assume that our problem ex- 
hibits separation of time scales; the amplitude changes 
on a typical time that is much longer than the period of 
oscillations. Hence, we can substitute the expressions for 
8 and Ku into the equation of motion, cf. Eq. (|33|l , and 
operate on both sides of the equation with the operator 

■§7 Jo"" {"} s i n ( wv t)dt to obtain 

- vw 2 A/2 + vw 4 A 3 /8 - vw 2 AA' 2 /8 + vA" /2 

- vw 2 A 2 A"/8 - 3vA' 2 A"/8 = (49) 
fwv 



y(x, t) = A(x + vt) sin[w(x + vt)} 



(45) 



2tt 

fwv 
~2^~ 



K}j{A(x + vt) sin[u>(a; + vt)]} sm(wvt)dt 



K H {A(x + vt) sm[w(x + vt)}} sxn(wvt)dt 
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This is a highly non-trivial integro-differential equation 
for the time evolution of the amplitude A. We expect 
that after a transient the amplitude saturates to a fixed 
value. Therefore, we set all the derivatives to zero and 
the amplitude equation reduces to 

- vw 2 A/2 + vw 4 A 3 /8 = (50) 



fwv 
~2^~ 

fwv 

~2^r 



Kjj {As'm[w(x + vt)]} am(wvt)dt 
V Kif {Asm[w(x + vt)]} sin(wvt)dt 



If we consider one of the control parameters slightly 
above its critical value, e.g. the velocity, and expand all 
the terms around the critical point we obtain 

- w 2 osc A/2 - e„(w 2 sc + 2v osc w osc d v w osc )A/2 

+ w 4 sc A 3 /8 = -fAK^{Asm[ Wosc (x + v osc t)}}/2v osc 

+ a(v osc )e v A + f3(v osc )A 3 . 

Here e v = v ~ Vosc is the critical parameter, a(v osc ) 
is the coefficient of a critical linear term related to 
d v K${A sm[w osc (x + v osc t)]} , (3(v osc ) is the coefficient 
of the non-critical cubic term, and we neglected terms of 
0(e v A 3 ). The linear non-critical terms cancel out since 
they are just the equation to O(A). Therefore, we are 
left with 







[(«>. 



2 

osc 
4 

osc 



2v, 



OSC^OSC 



+ [-w 4 osc /8 + f3( Vo 



cd v W 0: 



; )/2 + a{v osc )]e v A 



(51) 

where d v w osc > since as we increase the velocity we 
increase the energy flow to the crack tip which requires 
more crack surface created per unit time. The instability 
is tantamount to a positive linear coefficient. We assume 
that the cubic term leads to saturation. It follows that 



A = 



(^OSC 

+ 2^ osc w OS cG> 1 ,w osc )/2 + a{v osc ) 1/2 



(52) 



wt sc /8-P(v osc ) 

This result shows that the critical exponent with respect 
to e v is 1/2, in agreement with the experimental data0. 
A similar calculation can be presented for the dependence 
of A on eat = (AT — AT osc )/AT osc . Such a calculation 
result with the same exponent 1/2 as observed in the 
experiment. 

Two comments are in order. First, we derived the am- 
plitude equation from the dynamics of the tip, Eq. ({3*3"jl . 
rather than guess it as in previous works. More stamina 
can in principle lead to an actual calculation of the last 
term in this equation. Second, we have projected the full 
amplitude equation onto its asymmetric part. Projecting 
onto the symmetric (cosine) part yields the equation 

vwA' - vw 3 A 2 A'/2 - vwA' 3 /A - vwAA'A" /A = (53) 

— 

sin[u>(x + vt)]} cos(wvt)dt 

27r Jo 

2-7T 

- [ "™ Kf){A{x + vt) sxn[w(x + vt)]} cos(wvt)dt . 

2n Jn 



Discarding again, for the stationary state, all the deriva- 
tives, we see that we do not gain any new information 
about the stationary amplitude. On the other hand we 
learn that also the last term is a pure sine function, since 

it has to vanish at the critical point exactly like the Kj X j 
term. 



IV. CONCLUDING REMARKS 

The main point of departure of our theory from pre- 
vious ones is that we employ, in addition to the two- 
dimensional linear elasticity part, the dynamical crack- 
tip propagation law suggested in Ref . [9( . Using this dy- 
namical law we first derived a stability criterion for the 
straight crack propagation which is identical to a previ- 
ously suggested criterion 0, . We then extended the 
analysis to the evolution of the crack shape just above 
the onset of the oscillatory instability and showed that 
the dynamical equation has a stationary sinusoidal so- 
lution with a theoretically calculated wavenumber. We 
presented a quantitative comparison with the experimen- 
tal data for a temperature field that is characterized both 
by the spatial separation between the thermal baths h 
and the thermal diffusion length dt^. Our results agree 
rather well with the experiments p|. 

From the conceptual point of view we have offered a 
successful way to decompose the problem into a singular 
and a non-singular parts. This decomposition enabled us 
to derive an expression for Kn to leading order in the am- 
plitude of the oscillations that depended only on the fac- 
torization of one Wiener- Hopf kernel. This factorization 
is done by applying the method of Pade approximants 
suggested recently [l(j]- Finally, we showed how the dy- 
namical tip propagation law translates to an amplitude 
equation for the oscillatory solution. This equation re- 
sulted in calculated critical exponents of the transition, 
in agreement with the measured ones. The success of 
the dynamical theory based on the law of tip propaga- 
tion lends strong support to this law, at least in these 
quasi-static conditions. One should stress at this point 
that the analysis considered the temperature field as ef- 
fectively frozen. The oscillatory nature of the crack has 
very little to do with the temperature dynamics. This 
cannot be expected to remain valid for larger amplitudes 
of oscillations since the boundary condition restrict the 
temperature level sets to be normal to the crack. Thus at 
some point the dynamics of the termperature field must 
enter the discussion, potentially leading to new dynamic 
instabilities including chaos and disorder. 

In fact, the conclusion of this study appears to be that 
in the quasi-static conditions the assumption of small 
scale yielding holds, making it sufficient to solve the lin- 
ear elasticity problem, coupled to a correct law of motion 
that dictates how the tip propagates. It would be inter- 
esting to try to apply this or similar laws to other con- 
texts in which the quasi-static problem can be solved, 
but where the absence of an accepted propagation law 
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has led to a number of possible evolutions 0, 0] . We 
expect however that in truly dynamical crack propaga- 
tion new theoretical concepts need to be developed in 
order to reach a similar level of calculation of experimen- 
tal observations. 
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APPENDIX A: SLIGHTLY CURVED CRACKS 

The aim of this appendix is to derive Eq. (|39l) and 
Eq. 102). Assume that we have a mode I load cr yy (x, 0) 
on a semi-infinite crack whose tip is at x = 0. Let us first 
find the normal opening stress T n (x,y(x,t)) and tangen- 
tial shearing stress T t (x,y(x,t)) on any small deviation 
y(x,t) from straight crack in terms of a yy (x 1 0). The 
cartesian components of the stress tensor field are related 
to polar components according to 

&xx = &rr cos 2 8 + agg sin 2 9 — a T g sin 29 
G yy — a rr sin 2 8 + agg cos 2 9 + a r g sin 29 

<Jxy = (tfVr — CTgg) sin {29/2) + <J r g COS 28 . (Al) 

Here 8 is the local angle made by the tangent to the crack 
and the x-axis. These relations can be expanded to first 
order in 8 ~ y'(x, t) and then inverted to yield 

T n (x, y(x, t)) = <rgg(x, y(x, t)) (A2) 
= <Tyvfa v( x , *)) ~ 2 y'( x , t)(T X y[x, y(x, t)) 

T t (x, y(x, t)) = <r r g(x, y(x, t)) = a xy (x, y(x, t)) 
+y'(x, t) \a yy (x, y (x, t)) - a xx (x, y (x, t))] . 

Expanding the cartesian components to first order in 
y(x,t) 

cr yy (x,y(x,t)) = a yy (x, 0) + d y a yy (x, 0) y(x,t) 

(x,y{x,t)) 0) + d y a xy (x, 0) y(x, t) (A3) 

and using the relation 

d y a xy {x,0) = dy{-d x d v x) = ~d x d y d y x = -a xx (x,0) 

(A4) 

we end up with 

T n (x,y(x,t)) = a yy (x,0) 

T t (x,y(x,t)) = y'(x,t)((Xyy(x, 0) — <r xx (x,0)) 

- y(x,ty xx (x^) (A5) 



where we used the symmetry of the problem to set 

&xy 

0) = 0; a xx (x, 0) can be calculated from the knowl- 
edge of the boundary condition <J yy (x, 0). 

The problem we should solve now is formulated as fol- 
lows; Given the following crack configuration and load- 
ing conditions: (i) a semi- infinite crack whose shape is de- 
scribed by a small deviation y(x, t) from a straight crack 
configuration in an infinite two-dimensional domain (ii) 
a normal opening load T n {x,y{x,t)) and a shear load 
T t (x,y(x,t)) at any point on the crack, what are the 
mixed mode stress intensity factors? 

A version of this problem was treated completely by 
Cotterell and RicejlJ. They have found that the stress 
intensity factors for a finite slightly curved crack extend- 
ing from —a to a, where the deviation vanishes at both 
tips, are given by 

1 f a la + x 

Ki-iKu = \ dx [qi(x)-iqn(x)]\ (A6) 

V™ J -a \ a - x 

where qi (x) and qn (x) were derived explicitly in . In 
order to adapt this result to our semi-infinite configura- 
tion we should fix one tip of the crack to x = 0, take the 
limit where the other tip goes to — oo and finally shift the 
origin, 

y(x,t) = y(x,t) - y(0,t) , (A7) 

such that the deviation vanishes as the tip of the crack 
at any time. This adaptation yields 

Kl -iKn = f-f d ^Mf3iim (A8) 

V 7T J_ 00 sJ-X 

where 

qi = T n - ^'(0, t)T t + y(x, t)T[ + 2y{x, t)T t 



qn = T t + y{x,t)T^ + -y'{Q,t)T n (A9) 

Using the derived loading conditions of Eq. (|A5|I and 
applying these results to our ansatz for y(x,t) we obtain 

qi = <j yy {x,0) + 0(A 2 ) (A10) 



q n = {Asin[w(x + vt)](a yy (x,0) - a xx (x,Q))}' 

-A s'm(wvt)a' yy (x,0) + -Aw cos(wvt)a yy (x,0) + 0(A 3 ) . 

The last step is to relate a xx (x, 0) to the boundary 
condition a yy (x, 0). According to the complex potentials 
method |l3j we have for a semi-infinite straight crack 

o yy + a xx = 4Re[$(z)] (All) 

with 

= 1 f° dx[a yy (x,0) - i(Jxy(x,0)}y/=X 
2-k^/z J_ 00 z-x 

(A12) 
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Since in our case there is no shear loading on the straight 
crack we obtain 

qi = a yy (x,0) + 0{A 2 ) (A13) 
qu = {2Aain[w(x + vt)]a yy (x,0)}' — Asm(wvt)a' yy (x,0) 



+ — Aw cos(wvt) a yy (x, 0) +0(A 3 ) 



which becomes, upon substitution into Eq. (|A8|) . the 
RHS of Eq. (|42H and by setting t — and integrating by 
parts gives Eqs. 
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